
#include "StatefulRhoMaterial.h"

registerMooseObject("TrilobitaApp", StatefulRhoMaterial);

InputParameters
StatefulRhoMaterial::validParams()
{
    InputParameters params = Material::validParams();
    //add regular thermal properties
    params.addRequiredParam<Real>("rhov", "density of original material");
    params.addRequiredParam<Real>("rhoc", "density of fully charred material");
    //add special parameters to compute pyrolysis Arrhenius relation
    params.addRequiredParam<Real>("T_begin", "temperature at which the reaction begins");
    params.addRequiredParam<Real>("precofficient", "pre-exponential factor");
    params.addRequiredParam<Real>("ER", "E/R");
    params.addRequiredParam<Real>("reaction_level", "Reaction level of pyrolysis,normally 1,2 or 3");
    //add coupled variables
    params.addRequiredCoupledVar("temperature", "Coupled Temperature");
    return params;
}

Real StatefulRhoMaterial::Source(Real rho, Real T, Real T_begin, Real rhoc, Real rhov, Real precoff, Real ER, Real reaction_level)
{
    if (T > T_begin)
    {
        if ((rho - rhoc) / (rhov - rhoc) < 0.01)
        {
            return 0;
        }
        else
        {
            Real tmp = (rho - rhoc) / rhov;
            return precoff * exp(-ER / T) * rhov * pow(tmp, reaction_level);
        }
    }
    else
    {
        return 0;
    }
}

StatefulRhoMaterial::StatefulRhoMaterial(const InputParameters& parameters)
    : Material(parameters),
    _T_begin(declareProperty<Real>("T_begin")),
    _precofficient(declareProperty<Real>("precofficient")),
    _ER(declareProperty<Real>("_ER")),
    _reaction_level (declareProperty<Real>("reaction_level")),
    _rhov(declareProperty<Real>("rhov")),
    _rhoc(declareProperty<Real>("rhoc")),
    _rhos(declareProperty<Real>("rhos")),
    _pyrolysis_source(declareProperty<Real>("pyrolysis_source")),
    _rhos_old(getMaterialPropertyOld<Real>("rhos")),
    _T_value(coupledValue("temperature"))
{
    _rho_v_value = getParam<Real>("rhov");
    _rho_c_value = getParam<Real>("rhoc");
    _T_begin_value = getParam<Real>("T_begin");
    _precoff_value = getParam<Real>("precofficient");
    _ER_value = getParam<Real>("ER");
    _reaction_level_value = getParam<Real>("reaction_level");
    
}

void
StatefulRhoMaterial::initQpStatefulProperties()
{
    _rhov[_qp] = _rho_v_value;
    _rhoc[_qp] = _rho_c_value;
    _T_begin[_qp] = _T_begin_value;
    _precofficient[_qp] = _precoff_value;
    _ER[_qp] = _ER_value;
    _reaction_level[_qp] = _reaction_level_value;
    _rhos[_qp] = _rho_v_value;
}

void
StatefulRhoMaterial::computeQpProperties()
{
    _rhov[_qp] = _rho_v_value;
    _rhoc[_qp] = _rho_c_value;
    _T_begin[_qp] = _T_begin_value;
    _precofficient[_qp] = _precoff_value;
    _ER[_qp] = _ER_value;
    _reaction_level[_qp] = _reaction_level_value;
    Real source = Source(_rhos_old[_qp], _T_value[_qp], _T_begin_value, _rho_c_value, _rho_v_value, _precoff_value, _ER_value, _reaction_level_value);
    _rhos[_qp] = _rhos_old[_qp] - source * _dt;
       
}
